% This function conducts VAR analysis with long run restriction
% Blanchard-Quah (imposing the long-run effect of the second
% shock on the level of the first variable is 0, while the first 
% variable is a growth rate). Restricted to bivariate case.
% x: T*N data
% p: choice of lag
% J: how many periods of IRFs you want
% =====================================================
% By putting cons growth and cointegration between cons
% and corporate earning (1947Q1-2005Q4), p=5, J=80, I 
% replicate Hansen-Heaton-Li (2008 JPE?, Figure 3
% =====================================================
function [w,u,irf1,irf2,cum1,cum2]=HHL(z,p,J)
[T,N]=size(z);

% Reduced form VAR
y=z(p+1:end,:);
% regressors
x=z(p:end-1,:);
i=1;
while i<p
    x=[x,z(p-i:end-1-i,:)];
    i=i+1;
end
x=[x,ones(T-p,1)];
A0=(x'*x)\(x'*y); % VAR coefficients
u=y-x*A0; % residual
sigmau=u'*u/(T-p);

A=eye(N);
i=1;
while i<p+1
    A=A-(A0(N*(i-1)+1:N*i,1:N))';% A(1)
    i=i+1;
end

tempv=inv(A)*sigmau*inv(A');
theta=chol(tempv,'lower');
% Long-run restriction
b=A*theta; %B0^(-1)
w=(b*u')';

% impulse response
A_comp=(A0(1:N,1:N))';
i=2;
while i<p+1
    A_comp=[A_comp,(A0((i-1)*N+1:i*N,1:N))'];
    i=i+1;
end
A_comp=[A_comp;eye(N*(p-1)),zeros(N*(p-1),N)];
B_comp=[b,zeros(N,N*(p-1));zeros((p-1)*N,p*N)];

irf1=zeros(J,1);
irf2=zeros(J,1);
for i=1:J
    irf1_temp=A_comp^(i-1)*B_comp*[1;0;zeros(N*(p-1),1)];
    irf1(i,1)=irf1_temp(1);
    irf2_temp=A_comp^(i-1)*B_comp*[0;1;zeros(N*(p-1),1)];
    irf2(i,1)=irf2_temp(1);
end

cum1=zeros(J,1);
cum2=zeros(J,2);
for i=1:J
    cum1(i,1)=sum(irf1(1:i,1));
    cum2(i,1)=sum(irf2(1:i,1));
end
end